# -*- coding: utf-8 -*-
"""

PBAR_ExamineCalculateFeaturesZspecBasic.py

Test CalculateFeaturesZspecBasicInterp.py and CalculateFeaturesZspecBasic.py

Created on Thu May 01 15:27:24 2014

@author: jkwong
"""

import PBAR_Zspec, PBAR_Cargo
import os
import numpy as np
reload(PBAR_Zspec)
reload(PBAR_Cargo)

dataPath = r'E:\PBAR\data4\BasicScanCargo'
dataPathSW = r'E:\PBAR\data4\BasicScansStandardWidth'

datasetDescription = PBAR_Zspec.ReadCargoDataDescriptionFile(r'E:\PBAR\data4\CargoSet2.txt')

(goodZspecMask, badZspecMask, goodZspecNameList, badZspecNameList, \
        goodZspecIndices, badZspecIndices) = PBAR_Zspec.ZspecDetectorLists2ndSet()

energy = np.array([7.50E+00,1.05E+01,1.35E+01,1.65E+01,1.95E+01,2.25E+01,2.55E+01, \
    2.85E+01,3.15E+01,3.45E+01,3.75E+01,4.05E+01,4.35E+01,4.65E+01,4.95E+01,\
    5.25E+01,5.55E+01,5.85E+01,6.15E+01,6.45E+01,6.75E+01,7.05E+01,7.35E+01, \
    7.65E+01,7.95E+01,8.25E+01,8.55E+01,8.85E+01,9.15E+01,9.45E+01,9.75E+01,1.01E+02])

i = 0

filename = '%s-FDFC-All_SW.npy' %datasetDescription['scanID'][i]
fullfilename = os.path.join(dataPathSW, filename)
print('Loading %s' %fullfilename)

(energy, dat) =  PBAR_Zspec.ReadZspecBasicScanNumpy(fullfilename)
features = PBAR_Zspec.CalculateFeaturesZspecBasic(dat, energy, 7)
featuresInterp = PBAR_Zspec.CalculateFeaturesZspecBasicInterp(dat, energy, 7)


# manually calculate some values
ch = 50
t = 1000
binMeanTemp = np.sum(dat[t, ch,:]* energy) /np.float( dat[t,ch,:].sum())
# calculate interpreted value
energyFine = np.arange(energy[0], energy[-1], 1)
lowerBinThreshold = np.argmin( np.abs( energyFine - 7))
datFine = np.interp(energyFine, energy, dat[t,ch,:])
binMeanFineTemp = np.sum(datFine* energyFine) /np.float( datFine.sum())

print('Manually calculated value: %3.6f' %binMeanTemp)
print('Code calculated value: %3.6f' %features['binMean'][t, ch])
print('Manually calculated value, interp: %3.6f' %binMeanFineTemp)
print('Code calculated value, interp: %3.6f' %featuresInterp['binMean'][t, ch])

plt.figure()
plt.grid()
plt.plot(energy, dat[t, ch,:], 'o-b', alpha = 0.5, markersize = 10, linewidth = 4)
plt.plot(energyFine, datFine, 'xg', alpha = 0.5, markersize = 10)
plt.xlabel('Bin')
plt.ylabel('Count')

# plot the time profile for each feature
# - interpt and original
# Objective - show that the interpreted values are okay

featuresPlotList = features.keys()

featuresPlotList = ['binMean', 'binSTD', 'binSkew', 'binKur', 'multibin_20_ratio']

for featureName in featuresPlotList:
    plt.figure()
    plt.grid()
    for detIndex in [76]:
        plt.plot(features[featureName][:,detIndex], linestyle = '-', label = 'Normal, ch %d' %(detIndex+1), alpha = 0.5, linewidth = 2)
        plt.plot(featuresInterp[featureName][:,detIndex], color = 'r', linestyle = ':', label = 'Interp, ch %d' %(detIndex+1), alpha = 0.5, linewidth = 2)
    plt.legend()
    plt.title(featureName)

featureName = 'binKur'
plt.figure()
plt.imshow(features[featureName].T, interpolation = 'nearest', aspect='auto', cmap = cm.Greys_r)
plt.title('Feature: %s, Regular' %featureName)

plt.figure()
plt.imshow(featuresInterp[featureName].T, interpolation = 'nearest', aspect='auto', cmap = cm.Greys_r)
plt.title('Feature: %s, Interp' %featureName)
